An example of SSH reconstruction has been produced in the "example_data_oi.ipynb" notebook. Here, an example of data evaluation is proposed. The notebook is structured as follow: 1) reading of reference and reconstructed SSH fields, 2) make field on similar spatio-temporal grid and 3) comparison of reconstrusted and refernce SSH fields (statistical/spectral comparison)
import xarray as xr
import cftime
import geoviews as gv
import matplotlib.pylab as plt
from matplotlib.ticker import ScalarFormatter
import numpy as np
from datetime import datetime, timedelta
import numpy
import pyinterp
import holoviews as hv
import xrft
from dask.diagnostics import ProgressBar
gv.extension('bokeh')
from matplotlib import rcParams
rcParams['font.sans-serif'] = 'TeX Gyre Heros'
dc_ref = xr.open_mfdataset('./dc_ref/*.nc', combine='nested', concat_dim='time')
dc_ref
dc_reconstruction = xr.open_mfdataset('ssh_rec.nc', combine='by_coords')
dc_reconstruction
A regridding on the similar spato-temporal grid s required for the comparison. Here we have mde the choice to "temporally degrade" the reference SSH field into daily mean sample and interpolate the reconstructed field onto this new reference spatio-temporal grid. For this, the pyinterp package is used
dc_ref_sample = dc_ref.resample(time='1D').mean()
dc_ref_sample
x_axis = pyinterp.Axis(dc_reconstruction["lon"][:], is_circle=False)
y_axis = pyinterp.Axis(dc_reconstruction["lat"][:])
z_axis = pyinterp.TemporalAxis(dc_reconstruction["time"][:])
ssh_rec = dc_reconstruction["ssh_rec"][:].T
grid = pyinterp.Grid3D(x_axis, y_axis, z_axis, ssh_rec.data)
mx, my, mz = np.meshgrid(
dc_ref_sample.lon.values,
dc_ref_sample.lat.values,
z_axis.safe_cast(dc_ref_sample.time.values),
indexing="ij")
ssh_rec = pyinterp.trivariate(grid,
mx.flatten(),
my.flatten(),
mz.flatten(),
bounds_error=True).reshape(mx.shape).T
dc_reconstruction_interp = xr.Dataset({'sossheig' : (('time', 'lat', 'lon'), ssh_rec)},
coords={'time': dc_ref_sample.time.values,
'lon': dc_ref_sample.lon.values,
'lat': dc_ref_sample.lat.values,
})
dc_reconstruction_interp
time_selection = '2013-01-01T23:00:00'
plt.figure(figsize=(15, 5))
plt.subplot(121)
dc_reconstruction_interp.sossheig.sel(time=time_selection, method='nearest').plot(vmin=-0.2, vmax=1, cmap='gist_stern')
plt.title('Reconstruction')
plt.subplot(122)
dc_ref_sample.sossheig.sel(time=time_selection, method='nearest').plot(vmin=-0.2, vmax=1, cmap='gist_stern')
plt.title('Reference')
# plt.savefig('example_ssh.png')
# SSH reconstruction resampled, otherwise too heavy for display
dataset2 = gv.Dataset(dc_reconstruction_interp.coarsen({'lon':6, 'lat': 6, 'time':6}, boundary="trim").mean(), ['lon', 'lat', 'time'], 'sossheig')
images2 = dataset2.to(gv.Image).redim(sossheig=dict(range=(-0.2, 1.)))
# SSH reference resampled, otherwise too heavy for display
dataset3 = gv.Dataset(dc_ref_sample.coarsen({'lon':6, 'lat': 6, 'time':6}, boundary="trim").mean(), ['lon', 'lat', 'time'], 'sossheig')
images3 = dataset3.to(gv.Image).redim(sossheig=dict(range=(-0.2, 1.)))
# delta SSH
delta_ssh = (dc_reconstruction_interp - dc_ref_sample).coarsen({'lon':6, 'lat': 6, 'time':6}, boundary="trim").mean()
delta_ssh = delta_ssh.rename({'sossheig': 'delta_ssh'})
dataset1 = gv.Dataset(delta_ssh, ['lon', 'lat', 'time'], 'delta_ssh')
images1 = dataset1.to(gv.Image).redim(delta_ssh=dict(range=(-0.3, 0.3)))
layout = hv.Layout(images3.opts(width=400, height=300, cmap='gist_stern', colorbar=True, title='SSH reference', toolbar='above') + images2.opts(width=400, height=300, cmap='gist_stern', colorbar=True, title='SSH reconstruction', toolbar='above') + images1.opts(width=400, height=300, cmap='coolwarm', colorbar=True, title='SSH reconstruction - reference')).cols(2)
layout